set more off
program drop _all
program define PANELsample, rclass

        syntax, studies(integer) estperstudy(integer) totalobs(integer) alpha(real) ///
		theta(real) obs(integer)
        
		// Remove existing variables
        drop _all

		//We first create the matrix to store the results of each study
		set matsize 10000
		matrix A = J(`totalobs',1,.)
		matrix B = J(`totalobs',1,.)
		matrix C = J(`totalobs',1,.)
		matrix D = J(`totalobs',1,.)
		
		forvalues i = 1/`studies' {
		scalar lambdai = 0.5+30*runiform()
		scalar alphai = `alpha'+2*rnormal()
		forvalues j = 1/`estperstudy' {
        // STEP ONE: Create the data for each study and estimate an effect
		clear
		set obs 100
		generate x = rnormal()
		// Note that each "study" has a difference error variance, causing the estimate 
		// of the effect to be estimated with varying degrees of precision
		// This is for random effects
			scalar lambdaij = lambdai+`theta'*runiform()
			scalar alphaij = alphai+0.5*rnormal()
			generate e = lambdaij*rnormal()
			generate y = 1 + alphaij*x + e
		quietly regress y x 
		scalar coef = _b[x]
		scalar secoef = _se[x]
		scalar tcoef = coef/secoef
		scalar ID = `i'
		scalar obsno = (`i'-1)*`estperstudy'+`j'
		
		matrix A[obsno,1] = coef
		matrix B[obsno,1] = secoef
		matrix C[obsno,1] = tcoef
		matrix D[obsno,1] = ID
		}
}	
		matrix bob = A,B,C,D
		svmat bob
		rename bob1 effect 
		rename bob2 seeffect 
		rename bob3 teffect 
		rename bob4 ID
		generate pet = (1/seeffect)

		// First run this program once to get the pre-publication study sample data
		// To get post-publication study sample data, uncomment one of the two sections
		// below.
		generate dummy = 1
		//replace dummy = 0 if abs(teffect) < 2
		//replace dummy = 0 if effect < 0
		by ID, sort: egen select = mean(dummy)
		replace effect = cond(select<0.65,.,effect)
		replace seeffect = cond(select<0.65,.,seeffect)
		replace teffect = cond(select<0.65,.,teffect)
		replace pet = cond(select<0.65,.,pet)

		metareg effect , wsse(seeffect) mm
		return scalar I2 = e(I2)
		
		summ effect, detail
		return scalar effectMED = r(p50)
		return scalar effectMIN = r(min)		
		return scalar effectP5 = r(p5)
		return scalar effectP95 = r(p95)
		return scalar effectMAX = r(max)
		return scalar N = r(N)
		
		summ teffect, detail
		return scalar teffectMED = r(p50)
		return scalar teffectMIN = r(min)		
		return scalar teffectP5 = r(p5)
		return scalar teffectP95 = r(p95)
		return scalar teffectMAX = r(max)
		
		gen absteffect=abs(teffect)
		gen sig=(absteffect>=2)
		summ sig
		return scalar pctsig = r(mean)
		
end
		
		
		
